\documentclass[10pt,a4paper]{article} 

\usepackage{ctex} % 中文支持
\usepackage[top=2.5cm, bottom=2.5cm, left=2.5cm, right=2.5cm]{geometry} % 页边距
\usepackage{amsmath, amssymb} % 数学公式与符号
\usepackage{graphicx}
\usepackage{pythonhighlight}
\usepackage{url} 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{titling}
\setlength{\droptitle}{-2cm} % 标题上移

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 使表格美观
\usepackage{array}
\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
\setlength\extrarowheight{3pt}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%文档的题目、作者与日期
\author{五六七 }
\title{体检指标因子分析 }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

\maketitle

\begin{abstract}
对一些体检者的生化指标进行因子分析。
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\setcounter{tocdepth}{2}
%\renewcommand\contentsname{目录}
%
%\renewcommand {\baselinestretch} {1.3}\normalsize 
%\tableofcontents 
%\renewcommand {\baselinestretch} {1.0}\normalsize


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{问题描述}
有25名健康人参与了一次体检，测得7项生化检验指标如下表。请对该数据进行因子分析。

\begin{table}[ht!]\centering
\caption{体检数据} \vspace{0.2cm}
\begin{tabular}{|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|M{1.5cm}|}\hline 
$x_1$ & $x_2$ & $x_3$ & $x_4$ & $x_5$ & $x_6$ & $x_7$  \\ \hline 
3.76	&	3.66	&	0.54	&	5.28	&	9.77	&	13.74	&	4.78	\\ \hline 
8.59	&	4.99	&	1.34	&	10.02	&	7.5	&	10.16	&	2.13	\\ \hline 
6.22	&	6.14	&	4.52	&	9.84	&	2.17	&	2.73	&	1.09	\\ \hline 
7.57	&	7.28	&	7.07	&	12.66	&	1.79	&	2.1	&	0.82	\\ \hline 
9.03	&	7.08	&	2.59	&	11.76	&	4.54	&	6.22	&	1.28	\\ \hline 
5.51	&	3.98	&	1.3	&	6.92	&	5.33	&	7.3	&	2.4	\\ \hline 
3.27	&	0.62	&	0.44	&	3.36	&	7.63	&	8.84	&	8.39	\\ \hline 
8.74	&	7	&	3.31	&	11.68	&	3.53	&	4.76	&	1.12	\\ \hline 
9.64	&	9.49	&	1.03	&	13.57	&	13.13	&	18.52	&	2.35	\\ \hline 
9.73	&	1.33	&	1	&	9.87	&	9.87	&	11.06	&	3.7	\\ \hline 
8.59	&	2.98	&	1.17	&	9.17	&	7.85	&	9.91	&	2.62	\\ \hline 
7.12	&	5.49	&	3.68	&	9.72	&	2.64	&	3.43	&	1.19	\\ \hline 
4.69	&	3.01	&	2.17	&	5.98	&	2.76	&	3.55	&	2.01	\\ \hline 
5.51	&	1.34	&	1.27	&	5.81	&	4.57	&	5.38	&	3.43	\\ \hline 
1.66	&	1.61	&	1.57	&	2.8	&	1.78	&	2.09	&	3.72	\\ \hline 
5.9	&	5.76	&	1.55	&	8.84	&	5.4	&	7.5	&	1.97	\\ \hline 
9.84	&	9.27	&	1.51	&	13.6	&	9.02	&	12.67	&	1.75	\\ \hline 
8.39	&	4.92	&	2.54	&	10.05	&	3.96	&	5.24	&	1.43	\\ \hline 
4.94	&	4.38	&	1.03	&	6.68	&	6.49	&	9.06	&	2.81	\\ \hline 
7.23	&	2.3	&	1.77	&	7.79	&	4.39	&	5.37	&	2.27	\\ \hline 
9.46	&	7.31	&	1.04	&	12	&	11.58	&	16.18	&	2.42	\\ \hline 
9.55	&	5.35	&	4.25	&	11.74	&	2.77	&	3.51	&	1.05	\\ \hline 
4.94	&	4.52	&	4.5	&	8.07	&	1.79	&	2.1	&	1.29	\\ \hline 
8.21	&	3.08	&	2.42	&	9.1	&	3.75	&	4.66	&	1.72	\\ \hline 
9.41	&	6.44	&	5.11	&	12.5	&	2.45	&	3.1	&	0.91	\\ \hline 
\end{tabular}
\end{table}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{建立模型}
\subsection{对原始数据进行标准化处理}
将上述表格的每列数据都减去均值再除以标准差，计算公式为
\begin{eqnarray}
\tilde{x}_{i,j} = \frac{x_{i,j}-\mu_j}{s_j}, \,\,\, i=1,2,\cdots,25, \,\, j=1,2,\cdots,7. 
\end{eqnarray}
其中每列的均值和标准差的计算公式为 
\begin{eqnarray}
\mu_j &=& \frac{x_{1,j} + x_{2,j} + \cdots + x_{25,j}}{25}, \\ 
\sigma_j &=& \sqrt{\frac{1}{24}\left[ (x_{1,j}-\mu_j)^2 + (x_{2,j}-\mu_j)^2 + \cdots + (x_{25,j}-\mu_j)^2 \right] }. 
\end{eqnarray}

\subsection{计算任意两列之间的相关系数以及所有列的相关系数矩阵}
我们现在有7列数据，计算人意两列数据之间的相关系数的公式为
\begin{eqnarray}
r_{i,j} = \frac{1}{24} \left[ \tilde{x}_{1,i}\tilde{x}_{1,j} + \tilde{x}_{2,i}\tilde{x}_{2,j} + \cdots + \tilde{x}_{25,i}\tilde{x}_{25,j}  \right]. 
\end{eqnarray}
因此所有列的相关系数矩阵可以写为
\begin{eqnarray}
R=\frac{1}{24}\tilde{X}^t\tilde{X} = (r_{i,j})_{7\times 7},  
\end{eqnarray}
其中 $\tilde{X} = (\tilde{x}_{i,j})_{7\times 7}$ 是标准化处理后的数据矩阵。


\subsection{计算初等载荷矩阵}
相关系数矩阵是一个实对称矩阵，所以它的特征值都是实数，而且存在正交矩阵 $P$ 使得
所有列的相关系数矩阵正交相似于对角阵，并设已经从大到小排列， 
\begin{eqnarray}
P^{-1}RP = \begin{bmatrix} 
\lambda_1 & 0 & 0 & \cdots & 0 \\ 
0 & \lambda_2 & 0 & \cdots & 0 \\ 
0 & 0 & \lambda_0 & \cdots & 0 \\ 
\vdots & \vdots &\vdots & & \vdots \\ 
0 & 0 & 0 & \cdots & \lambda_7 \\ 
\end{bmatrix}. 
\end{eqnarray}
因为矩阵 $R$ 是半正定的，所以特征值都是非负的，即有 
\begin{eqnarray}
\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge \cdots \ge \lambda_7 \ge 0. 
\end{eqnarray}

记正交矩阵 $P$ 的列向量分别为 
\begin{eqnarray}
\vec{u}_1 = \begin{bmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{7,1} \end{bmatrix}, \,\, 
\vec{u}_2 = \begin{bmatrix} u_{1,2} \\ u_{2,2} \\ \vdots \\ u_{7,2} \end{bmatrix}, \,\, 
\cdots, \,\, 
\vec{u}_7 = \begin{bmatrix} u_{1,7} \\ u_{2,7} \\ \vdots \\ u_{7,7} \end{bmatrix}. 
\end{eqnarray}
使用主成分分析法，得到初等载荷矩阵
\begin{eqnarray}
A_1 = 
\left[ 
\sqrt{\lambda_1} \vec{u}_1, 
\sqrt{\lambda_2} \vec{u}_2, 
\cdots, \,\, 
\sqrt{\lambda_7} \vec{u}_7 
\right]. 
\end{eqnarray}
则有 $R = A_1A_1^t$. 


\subsection{构造因子模型}

计算特征值的贡献率和累计贡献率，确定因子的个数，然后构造因子模型 
\begin{eqnarray}
\vec{x} = \vec{\mu} + A\vec{f} + \vec{\varepsilon},
\end{eqnarray}
其中 $\vec{\mu}$ 是数据的期望向量，$\vec{f}$ 是公共因子向量，$\vec{\varepsilon}$ 是特殊因子向量。
矩阵 $A=(a_{ij})$ 是因子载荷矩阵，$a_{ij}$ 是变量 $x_i$ 在 公共因子 $f_j$ 上的载荷，反映了因子 $f_j$ 对 $x_i$ 的重要度。模型的假设有以下几条：
\begin{enumerate}
\item  公共因子 $f_j$ 互不相关且具有单位方差，
\item  特殊因子 $\varepsilon_i$ 互不相关且与 $f_j$ 互不相关，协方差矩阵 $\text{cov}(\vec{\varepsilon},
\vec{\varepsilon})=\Psi$ 为对角阵。
\end{enumerate}

在这些假设下，计算协方差矩阵可得
\begin{eqnarray}
\text{cov}(\vec{x},\vec{x}) &=& AA^t + \Psi, \\ 
\text{cov}(\vec{x},\vec{f}) &=& A. 
\end{eqnarray}

\subsection{本题中的因子模型}

在本题中，首先计算相关系数矩阵的特征值和特征向量。特征值从大到小排列如下，
\begin{table}[ht!]\centering
\caption{特征值以及各因子的贡献}\vspace{0.2cm}
\begin{tabular}{|c|c|c|c|c|c|c|c|} \hline 
特征值 &3.3952& 2.8063& 0.4365& 0.2762& 0.0812& 0.0042& 0.0004 \\ \hline 
贡献率 &0.485& 0.4009& 0.0624& 0.0395& 0.0116& 0.0006& 0.0001 \\ \hline 
累计贡献率 &0.485& 0.8859& {\color{red}0.9483}& 0.9877& 0.9993& 0.9999& 1.0000 \\ \hline 
\end{tabular}
\end{table}

因此选取 $m=3$ 个主因子，可以使得累计贡献率为 $0.95$. 因子模型为
\begin{eqnarray}
\begin{bmatrix} \tilde{x}_1 \\  \tilde{x}_2 \\ \vdots \\  \tilde{x}_7 \end{bmatrix}
= \begin{bmatrix} 
\alpha_{11} & \alpha_{12} & \alpha_{13} \\
\alpha_{21} & \alpha_{22} & \alpha_{23} \\
\vdots & \vdots & \vdots \\ 
\alpha_{71} & \alpha_{72} & \alpha_{73} \\
 \end{bmatrix} 
\begin{bmatrix} f_1 \\ f_2 \\ f_3  \end{bmatrix}  
+ 
\begin{bmatrix} \varepsilon_1 \\  \varepsilon_2 \\ \vdots \\  \varepsilon_7 \end{bmatrix}, 
\end{eqnarray}
其中因为数据 $\tilde{x}$ 已经标准化，所以期望向量 $\vec{\mu}$ 为零。
因子载荷矩阵为 
\begin{eqnarray}
A_1 = \begin{bmatrix} 
\alpha_{11} & \alpha_{12} & \alpha_{13} \\
\alpha_{21} & \alpha_{22} & \alpha_{23} \\
\vdots & \vdots & \vdots \\ 
\alpha_{71} & \alpha_{72} & \alpha_{73} \\
 \end{bmatrix} 
= \begin{bmatrix} 
\sqrt{\lambda_1}u_{1,1} & \sqrt{\lambda_2}u_{1,2} & \sqrt{\lambda_3} u_{1,3} \\
\sqrt{\lambda_1}u_{2,1} & \sqrt{\lambda_2}u_{2,2} & \sqrt{\lambda_3} u_{2,3} \\ 
\vdots & \vdots & \vdots \\ 
\sqrt{\lambda_1}u_{7,1} & \sqrt{\lambda_2}u_{7,2} & \sqrt{\lambda_3} u_{7,3} \\ 
\end{bmatrix}
= \begin{bmatrix} 
-0.7465 &  -0.4893 &  -0.4426 \\ 
-0.7964 &  -0.3722 &   0.4598 \\ 
-0.7089 &   0.5973 &   0.0995 \\ 
-0.9105 &  -0.3887 &  -0.074  \\ 
0.2342 &  -0.9635 &   0.0193 \\ 
0.1771 &  -0.9717 &   0.1148 \\ 
0.8864 &  -0.2192 &   0.016  \\ 
\end{bmatrix}. 
\end{eqnarray}

\subsection{因子旋转}

为了使得因子载荷矩阵的结构简化，使其每列或每行的平方值趋于0和1两极，使用正交矩阵进行旋转。
调用 \verb+factor_analyzer+ 的库函数，可得旋转后的因子载荷矩阵为 
\begin{eqnarray}
A_2 = \begin{bmatrix} 
0.127  & {\color{red}0.9635}  &  0.2267 \\ 
0.0572  &  0.3772  &  {\color{red}0.9226} \\ 
-0.771  &   0.2406  &  0.3452 \\ 
-0.0076  &  0.8088  &  0.5867 \\ 
{\color{red}0.9873}  &  0.1527  &  0.0533 \\ 
{\color{red}0.9785}  &  0.1312  &  0.1655 \\ 
0.5113  & -0.498  &  -0.473 \\ 
\end{bmatrix}. 
\end{eqnarray}
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{编程计算}

载入数值计算包numpy, 数据预处理函数 zscore.  
\begin{python}
import numpy as np
from scipy.stats import zscore
\end{python}

载入数据，将数据标准化，求相关系数矩阵。
\begin{python}
c=np.loadtxt('ti12_2.txt')
n=c.shape[0]
d=zscore(c,ddof=1)
r=np.corrcoef(d.T)
\end{python}

计算相关系数矩阵的特征值和特征向量，求特征值的累加和。
\begin{python}
val,vec=np.linalg.eig(r)
cs=np.cumsum(val)
\end{python}

求每个特征值的贡献率，求累计贡献率。
\begin{python}
rate=val/cs[-1]
Trate=np.cumsum(rate)
\end{python}

打印特征值，打印贡献率，打印累计贡献率。
\begin{python}
print(val)
print(rate)
print(np.round(Trate,4))
\end{python}

选取前三列，得到旋转前的因子载荷矩阵。使用方差最大化方法，手工计算旋转矩阵，使得矩阵A1的列向量的分量的方差最大。打印旋转后的因子载荷矩阵。这里还没完成。
\begin{python}
A0=vec*np.sqrt(val)
A=A0[:,0:3]
P=
A1=A@P
print(A1)
\end{python}

直接调用 \verb+factor_analyzer+ 包，构建因子分析模型，求解方差最大的因子载荷矩阵。提取载荷矩阵并打印。
\begin{python}
from factor_analyzer import FactorAnalyzer as FA
fa = FA(3,rotation='varimax')  
fa.fit(d)
A2=fa.loadings_
print(A2)
\end{python}

因子载荷矩阵的每列元素求平方和，得到信息贡献。每行元素求平方和，得到共同度。计算因子得分。
\begin{python}
gx=np.sum(A2**2, axis=0)
s2=np.sum(A2**2, axis=1)
ss=np.linalg.inv(np.diag(1-s2))
f=ss@A2@np.linalg.inv(A2.T@ss@A2)
print(np.round(s2,4))
print(np.round(gx,4))
\end{python}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{回答问题}
从旋转后的因子载荷矩阵可以看出，第一个因子是 $x_5$ 和 $x_6$ 的综合因子，第二个因子是 $x_1$, 第三个因子是 $x_2$. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{参考文献 }
\begin{thebibliography}{99}

%\bibitem{dingtongren} 丁同仁、李承治，常微分方程教程，高等教育出版社，2022年3月第三版。
\bibitem{sishoukui-2} 司守奎,孙玺菁. \emph{Python数学建模算法与应用}, 国防工业出版社. 2022年1月第1版. 
\bibitem{hexiaoqun-ara} 何晓群. \emph{应用回归分析(R语言版)}. 电子工业出版社. 2017年7月第1版. 
%\bibitem{dalgaard} Peter Dalgaard 著, 郝智恒等译. \emph{R语言统计入门}. 人民邮电出版社. 2014年6月第1版. 


\end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





